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1.  Introduction 


Non-crystalline  (amorphous)  ceramics  or  ceramic  glasses  are  used  in  a  variety  of  vital  Army 
personnel,  ground,  and  air  vehicle  applications  that  require  transparent  armor — it  is  ubiquitous  in 
tactical  vehicular  windshields  and  side  windows.  Ceramic  glass  is  inexpensive  and  is  formable 
into  large,  flat  plate  and  curved  shapes.  For  many  years  it  has  been  known  that  the  properties  of 
glass  can  be  modified  and  enhanced  through  compositional  modification,  chemical  strengthening, 
annealing,  and  process  control  of  melt  cooling.  Glass  ceramics,  the  controlled  crystallization  of 
nano-sized  single  crystals  in  a  glass  matrix,  offer  another  avenue  for  designed  and  enhanced 
property  modifications  for  transparent  and  opaque  material  applications.  In  addition,  certain 
glass  formulations  have  been  shown  to  exhibit  enhanced  performance  against  shaped-charge  jets 
(SCJs)  (figure  1)  and  other  ballistic  threats,  but  it  is  not  understood  why.  This  is  in  part  due  to 
various  short-  and  longer-range  atomic  structural  characteristics  including  atomic  free  volumes, 
cation  coordinations,  bridging  and  non-bridging  oxygen  (O)  atoms,  bonding  energies,  and 
nanoscale  order  characteristics  (short  and  longer  range)  that  are  difficult  or  impossible  to  quantify 
experimentally  for  ceramic  glass. 

In  contrast,  crystalline  ceramics  like  silicon  carbide  ( SiC ),  aluminum  oxynitride  ( AlON ),  and 
others  have  easily  characterizable  microstructures/mesostructures,  which  consist  of  assemblies  of 
individual  single-crystal  grains.  Ceramic  glasses,  on  the  other  hand,  do  not  have  a  conventional 
micro-  or  mesostructure,  as  it  is  understood  for  crystalline  ceramics.  However,  there  are 
microstructural-scale  variations  in  ceramic  glass  that  may  include  density  variations  from  atomic 
free  volume  variations  or  microporosity,  size  of  local  atomic  order,  defects  (inclusions,  large 
pores,  etc.),  and  others  yet  to  be  determined. 

The  interaction  of  a  stress/shock  wave  from  a  dynamic  impact  involves  many  structural  changes 
not  easily  characterized  by  conventional  equations  of  state  and  can  involve  reversible  and 
irreversible  densification  and  changes  in  bulk  short-range  order  structures  comparable  to  phase 
changes  in  crystalline  ceramics.  For  example,  in  simple  Hertzian  indentation  testing,  a  wide 
range  of  plastic  or  inelastic  deformation  mechanisms  have  been  observed  in  a  variety  of  glasses. 
Multiscale  computational  design  methodologies  (figure  2),  for  this  class  of  materials  will, 
nevertheless,  require  quantitative  and  possibly  statistically  based  descriptions  of  the  mesoscale, 
although  current  efforts  to  develop  such  models  have  fallen  far  short  of  this  goal. 
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(a) 


(b) 


Figure  1.  Enhanced  performance  of  SCJs  into  glass  (a)  test  configuration  for  glass  targets, 
and  (b)  penetration  vs.  time  for  several  targets,  after  Moran  et  al.  (7). 
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Figure  2.  A  multiscale  model  for  non-crystalline  ceramics  (glass). 
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Our  specific  long-term  research  goals  are  threefold: 


1.  Develop  molecular  dynamics  (MD)  process  models  for  a  series  of  chemically  substituted 
amorphous  silica  (a -SiOo  or  fused  silica)  materials  for  the  prediction  of  glass  elastic  properties 
assuming  completely  uniform  glass  “mesostructures.”  If  successful,  such  models  will  enable 
ab  initio  prediction  of  structure-property  relations  in  glass  that  will  be  validated  with 
experimentally  determined  elastic  properties. 

2.  Extend  the  MD  models  to  study  densification  of  chemically  substituted  a -Si02  materials 
under  high  pressures  (~60  GPa  Materials  in  Extreme  Dynamic  Environments  [MEDE]) 
relevant  to  ballistic  events  where  reversible  and  irreversible  density  changes  and  structural 
transformations  have  been  observed.  If  successful,  such  models  will  enable  ab  initio 
prediction  of  a -S1O2  compressibility,  kinetics,  and  “glass”  phase  transformations  that  will  be 
used  to  develop  equations  of  state  for  a-SiC>2  materials,  and  thus  form  a  direct  link  to  the 
continuum  scale. 

3.  Develop  a  fully  validated  multiscale  finite  element  computational  model  and  code 
incorporating  the  effects  of  reversible  and  irreversible  densification,  inelastic  deformation,  and 
overlain  by  a  spatiotemporally  evolving  population  of  growing  defects,  which  coalesce  and 
lead  to  ultimate  fracture  and  fragmentation.  It  is  envisioned  that  at  some  time  in  the  not  too 
distant  future,  fully  concurrent  multiscale  computational  finite  element  codes  will  be  used  by 
analysts  on  a  regular  basis  for  optimal  material  design. 

The  remainder  of  the  report  is  organized  as  follows.  General  program  objectives  are  outlined  in 
section  2,  and  the  approach  for  modeling  the  multiscale  behavior  of  glass  appears  in  section  3. 
Experimental  work  on  various  glasses  is  described  in  section  4,  which  also  includes  an  overview 
of  ballistic  experiments  and  fragmentation  studies  conducted  at  the  Ernst-Mach  Institute; 
section  4  describes  nanoindentation  tests  which  are  useful  for  computational  validation  of  the  MD 
(section  5)  and  peridynamics  (section  7)  modeling  efforts.  A  shock  Hugoniot  for  fused  silica  is 
determined  using  force-matching  techniques  and  described  in  section  5.  First  principles  quantum 
mechanical  methods  are  used  to  model  densification  and  bulk  modulus  variation  with  pressure  in 
section  6.  Since  inelastic  deformation  and  a  polyamorphic  transformation  occurs  simultaneously 
in  fused  silica,  a  model  is  developed  to  account  for  this  coupled  behavior  in  section  8  and 
compared  to  plate  impact  experiments.  A  short-term  conceptual  project  to  determine  an  effective 
experimental  and  theoretical  approach  to  model  and  characterize  the  role  of  glassy  materials  in 
resisting  ballistic  impact  was  conducted  by  Professor  Richard  Lehman,  Rutgers  University,  in 
section  9.  Section  10  summarizes  the  conclusions  of  this  progress  report. 
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2.  Program  Objectives 


The  long-term  research  goal  of  the  program  is  to  develop  a  concurrent  multiscale  computational 
finite  element  code  for  optimizing  or  enhancing  the  performance  of  various  glasses  against  SCJs; 
the  initial  work  focuses  on  pure  fused-silica  (a -Si02)  and  chemically  varied  a-Si02  materials. 

As  such,  this  objective  falls  squarely  within  the  purview  of  the  Weapons  and  Materials  Research 
Directorate  (WMRD),  since  multiscale  models  are  constitutive  models  (specific  to  a  particular 
material)  wherein  time-evolving  microstructural  changes,  such  as  microcrack  growth,  are  fully 
coupled  to  the  macroscale,  a  phenomenon  that  cannot  be  modeled  or  accounted  for  using  classical 
homogenization  methods.  A  more  immediate  research  objective  is  to  understand  why  certain 
chemically  substituted  a -Si02  materials  exhibit  enhanced  performance  in  the  defeat  of  SCJ  and 
other  ballistic  threats. 

Our  program  objectives  are  threefold: 

1.  Develop  MD  process  models  for  a  series  of  chemically  substituted  a-Si02  materials  for  the 
prediction  of  glass  elastic  properties.  This  glass  plays  an  important  role  in  many  technological 
applications  and  its  structure  has  been  inferred  from  neutron-diffraction,  nuclear  magnetic 
resonance,  and  small-angle  x-ray  scattering  (SAXS)  analysis  to  reveal  a  three-dimensional 
network  consisting  of  tetrahedrally  coordinated  silicon  (Si)  whose  structure  is  constant 
throughout  the  glass  and  defines  its  short-range  order  (SRO).  Long-range  disorder  in  the 
structure  is  manifested  by  a  seemingly  random  variation  in  the  Si-O-Si  bond  angle  in  adjacent 
tetrahedra.  Despite  the  intense  study  of  a -Si02  glass  over  the  last  several  decades,  much 
controversy  still  exists  on  the  best  method  to  model  (i.e.,  via  density  functional  theory,  MD, 
Monte  Carlo  methods,  or  master  equation  techniques)  this  archetypal  material  for  prediction  of 
elastic  properties,  diffusivity,  surface  interactions,  bond  angle  distribution,  polyamorphism, 
and  melt  solidification.  Current  models  that  appear  in  the  literature  are  often  not  fully 
validated  and  progress  towards  this  goal  will  be  made  when  model  predictions  of  elastic 
constants  for  a  series  of  chemically  substituted  a -Si02  glasses  agree  with  experimentally 
determined  constants. 

2.  Extend  the  MD  models  to  study  densification  of  the  chemically  substituted  a -Si02  materials 
under  high  pressures.  Since  long-range  order  in  glass  is  non-existent,  variations  in  the  SRO, 
and  intermediate  range  order  (IRO)  must  be  responsible  for  the  enhanced  performance 
observed  in  ballistic  tests  on  certain  a -Si02  glasses.  If  this  is  the  case,  it  may  be  possible  to 
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use  MD  models  to  predict  macroscopic  ballistic  performance.  Since  glass  is  subjected  to 
extreme  pressure  and  temperature  during  an  SCJ  event,  it  will  be  necessary  to  study  the 
relationship  between  compressibility,  kinetics,  and  phase  transitions  during  high-pressure 
densification  of  a -SiOo  glasses  as  manifested  by  changes  in  coordination  number,  ring  size, 
and  free  volume.  Progress  towards  this  goal  will  be  made  when  MD-derived  equations  of 
state  (EOSs)  agree  with  those  obtained  experimentally  via  diamond  anvil  press  and  plate 
impact  experiments. 

3.  Develop  a  fully  validated  multiscale  finite  element  computational  model  and  code  that 

incorporates  the  effects  of  reversible  and  irreversible  densification,  and  inelastic  deformation, 
overlain  with  a  spatiotemporally  evolving  population  of  defects  that  grow,  coalesce,  and  lead 
to  ultimate  fragmentation.  This  objective  will  develop  a  computational  framework  to  combine 
the  objectives  from  (1)  and  (2),  and  incorporate  the  influence  of  fracture  initiation,  growth, 
coalescence,  and  fragmentation  of  surface  and  volume  defects  in  glass  into  a  comprehensive 
concurrent  multiscale  finite  element  model  and  code.  Microcrack  initiation,  growth,  and 
coalescence  (sometimes  referred  to  as  failure  waves)  is  a  multiscale  phenomenon  that  bridges 
all  scales  in  a-SiC>2  glasses  despite  the  apparent  absence  of  a  structural  mesoscale  for  this 
class  of  materials  (figure  2).  Algorithms  for  the  development  of  fully  two-way  coupled 
multiscale  codes  are  in  their  infancy,  and  progress  on  this  objective  will  be  realized  with 
successful  development  and  implementation  of  a  consistent  scheme  for  coarse-graining 
localization  phenomena  such  as  fracture  failure  observed  in  glass. 


3.  Planned  Approach 


The  planned  approach  consists  of  three  components,  which  are  outlined  in  figure  3: 

1.  Validate  the  MD  models  for  a  series  of  chemically  substituted  a-SiOj  materials  for  the 
prediction  of  glass  elastic  properties.  Although  there  is  no  effort  within  WMRD  to  predict 
a-SiC>2  elastic  properties,  a  hierarchical  multiscale  modeling  effort  is  currently  underway 
within  WMRD,  which  is  focused  on  the  study  of  polycrystalline  (~200  pm  grain  size)  AlON 
and  validation  of  quantum  and  MD  predictions  of  anisotropic  elastic  constants  using  diamond 
anvil  cell  (DAC)  and  focused-ion-beam  (FIB)/scanning  electron  microscopy  (SEM) 
compression  tests  on  oriented  AlON  single  crystals  (14).  We  plan  to  use  MD  methods  (with 
possible  MD  coarse-graining)  to  simulate  glass  process  modeling  during  melt  solidification  by 
quenching  a  high-density,  high  temperature,  and  pressure  (with  possible  polyamorphic  phases) 
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melt  for  a  series  of  chemically  substituted  a -Si02  materials.  Next,  the  resulting 
room-temperature,  chemically  modified  structures  will  be  reversibly  deformed  to  predict  the 
elastic  properties  that  will  be  validated  with  experimentally  determined  elastic  properties. 


FYll 

Glass  (DSI) 


FY12 


I.  Glass  Fabrication 


Fabricate  pure  &  chemically 
'  substituted  o-S/02  mat1  Is  (fused 
silica) 


FY13 


FY14 


FY15 


II.  Glass  Characterization 

T 

•  Evaluate  a-Si02  samples  in 
ballistic  &SCJ  tests 

•  Gla  s  s  Tech  n  ol  ogy  Sh  ort  Cou  rse 
(October  2010) 

•  Experimental  determination  of 

•  Quantitative  description  of 
r — ^initial  defect  distributions  in  o- 

y  5/02  materials 

•  Radial  distribution  function 

•  Novel  metrics  to  predict 
ballistic  performance  based 
on  MD  model  results 

microstructural  scale  variations  in 
glass  (SR0  &  LR0) 

measurements 

III.  Quantum  &  Molecular 
Dynamics  of  Glass 


MD  process  models  of  pure  & 
chemically  substituted  a-Si02 
mat'ls  validated  by  SR0& 
bond  length  distributions  from 


V*  MD  process  models  for  densification 
of  pure  &  chemically  substituted  a- 
Si02  materials  for 
Analytical  EOS  development 
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IV.  Meso/Macroscale  Modeling 
and  Validation 


Validate  MD  model  predictions  of 
elastic  properties  w/  experiments 
Begin  evaluation  of  coarse-graining 
algorithms  (CZ->XFEM) 


*  Develop  analytical  EOS  based 
for  chemically  substituted  o- 

VSi02  materials 

•Validate  MD  model 
densification  predictions 
w/plate  impact  experiments 


•  Implement  EOS  into 
continuum  codes  &  validate 
w/analysis  &  experiments 


Concurrent  multiscale 
computational  code  for  o- 
Si02  materials 


V.  Fundamental  Aspects 
of  Glass  Deformation 
&  Failure 

V 

•  Coherent  gradient  sensing/high-speed 
r — 7  imaging  experiments  on  pure  &  r— 7 

y  chemically  substituted  o-S/02  matt's  y 

•Plate  impact  experiments  on  pure& 
chemically  substituted  o-S/02  mat'ls 

•  Fragmentation  experiments 
on  a-SiO 2  mat'ls 

▼ 

Implement  &  validate 
coarse-graining  algorithms^ 
for  mat1 1  failure  of  o-S/02  , 
mat'ls 

Figure  3.  Five-year  roadmap  consistent  with  the  WMRD  brittle  materials  program. 


2.  Validate  the  MD  models  for  densification  of  chemically  substituted  a-SiO 2  materials  under 
high  pressures.  Although  there  is  currently  no  effort  within  WMRD  to  predict  the  EOS  of 
chemically  substituted  a -Si02  materials,  MD  methods  have  been  used  to  predict  high- 
pressure  densification  in  these  materials.  MD  simulations  of  pure  a -Si02  materials  reveal  a 
Hugoniot  elastic  limit  (HEL)  of  about  10  GPa,  and  an  anomalous  maximum  in  compressibility 
at  around  3  GPa.  Experiments  where  samples  have  been  compressed  to  pressures  lower  than 
10  GPa  are  indistinguishable  from  the  original  material,  whereas  above  10  GPa,  materials  can 
sustain  an  irreversible  density  increase  from  10-20%  higher  than  the  starting  material, 
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although  there  is  controversy  as  to  whether  the  mechanism  is  due  to  irreversible  coordination 
defects  or  permanent  ring  size  modification.  In  contrast  to  the  behavior  of  a-SiC>2,  crystalline 
quartz  (a-SiOf),  undergoes  very  well-known  high  pressure  polymorphic  phase  transitions  into 
a  Coesite  phase  and  a  Stishovite  phase  (figure  4),  which  involve  changes  in  coordination  of  the 
Si  cation  from  four  to  six  0  atoms.  A  combination  of  diamond  anvil  press  and  plate  impact 
experiments  will  be  conducted  on  a  series  of  chemically  substituted  a-SiOo  materials  and 
compared  with  MD  densification  simulations  in  glass  in  order  to  understand  the  influence  of 
glass  modifiers  on  changes  in  the  shock  response  of  these  materials.  EOSs  for  a  subset  of 
promising  chemically  substituted  materials  will  be  developed  and  implemented  into  a 
continuum  code  to  determine  if  any  of  the  chemically  substituted  materials  exhibit  enhanced 
ballistic  performance. 


T  (°C) 


Figure  4.  Crystal  structures  of  quartz,  after  Frye  (2). 


3.  Develop  a  fully  validated  multiscale  finite  element  computational  model  and  code  that 

incorporates  the  effects  of  reversible  and  irreversible  densification,  and  inelastic  deformation, 
overlain  with  an  spatiotemporally  varying  population  of  defects  which  grow,  coalesce,  and 
lead  to  ultimate  material  fragmentation.  The  ultimate  research  objective  is  to  develop  a 
physics-based  multiscale  computational  finite  element  code  for  studying  densification,  and 
dynamic  fracture  in  non-crystalline  ceramics  (see  figure  2).  Atomistic  behavior  will  be  linked 
to  macroscopic  elastic  properties  and  densification  behavior  through  development  of  an  EOS 
from  first  principles  as  outlined  in  components  (1)  and  (2)  above.  At  this  stage,  what  remains 
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to  be  accomplished,  is  to  successfully  link,  in  a  concurrent  fashion,  multiscale  failure 
phenomena  in  a-Si02  materials  by  incorporating  the  important  role  that  pre-existing  surface 
and  volume  defects  have  on  the  microcrack  growth,  coalescence,  and  fragmentation  in  this 
class  of  materials.  Over  the  past  five  years,  the  first  author  has  also  been  directly  involved  in 
development  of  a  parallel,  concurrent  multiscale  code  for  heterogeneous  viscoelastic 
composites  (15)  under  the  auspices  of  a  U.S.  Army  Research  Laboratory  (ARL)/University  of 
Nebraska  cooperative  agreement,  which  will  be  leveraged  and  used  as  the  framework  for  the 
development  of  a  concurrent  multiscale  model  of  a -Si02  materials. 

The  chief  challenge  for  brittle  materials  is  to  correctly  account  for  the  growth  kinetics  of 
microcracks  in  a  multiscale  computational  environment.  The  propagation  of  free  internal 
boundaries  at  lower  scales  will  be  “coarse-grained”  to  higher  scales,  where  global  fracture 
failure  and  fragmentation  is  observed.  As  such,  coarse-graining  algorithms  will  need  to  be 
validated  through  continuum-scale  experiments  on  a -Si02  materials  that  measure  dynamic 
crack  propagation  speeds,  mixed-mode  failure,  and  crack  bifurcation  phenomena  using 
coherent  gradient  sensing  and  high-speed  imaging  techniques;  ARL  possesses  capabilities  for 
conducting  such  dynamic  fracture  experiments  in  a-Si02  materials  through  the  recent 
establishment  of  a  coherent  gradient  sensing/imaging  facility  funded  by  the  ongoing 
multiscale  modeling  effort  of  AlON .  Models  and  validation  of  the  initiation  and  propagation 
of  discrete  fractures  in  a-Si02  materials  should  transition  naturally  into  models  of 
fragmentation  and  comminution  for  behind-armor-debris  applications.  Fragmentation 
experiments  have  classically  been  conducted  using  dynamically  expanding  ring  experiments 
for  defining  fragment  size  versus  strain  rate  and  will  be  used  to  validate  computational  models 
of  fragmentation.  The  development  of  consistent  coarse-graining  algorithms  for  fracture  in 
materials,  which  is  associated  with  failure  and  loss  of  material  stability,  is  largely  unexplored 
and  is  the  primary  high-risk  goal  of  this  section. 
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4.  Experimental  Work  and  Background 


4.1  Background 

4.1.1  Compositions 

Silicate-based  ceramic  glasses  are  based  on  chemical  substitutions  into  a  Si()\  tetrahedral-based 
polymeric-like  structure;  fused  silica  is  an  amorphous  (non-crystalline)  form  of  pure  Si02.  Table 
1  lists  the  compositions  and  properties  of  typical  glasses. 


Table  1 .  Glass  compositions  in  %  and  selected  properties. 


Glass 

Si02 

Al2o3 

CaO 

B2o3 

Na20 

k2o 

p{g/ cm3) 

E  (GPa) 

V 

Borofloat 

80.5 

2.5 

0.02 

12.7 

3.55 

.64 

2.23 

62.3 

.207 

Starphire 

73.2 

1.44 

10.27 

- 

14.72 

.01 

2.51 

72.3 

.23 

Fused  Silica 

100 

- 

- 

- 

- 

- 

2.2 

73.0 

.17 

4.1.2  Structural  Characteristics  of  Ceramic  Glasses 

Simplistically,  the  predominant  macro-characteristics  (micron  and  larger)  can  be  a  variety  of 
defects  including  inclusions,  bubbles,  large  pores,  and  residual  compressive  or  tensile  stress.  The 
notion  of  an  array  of  crystalline  grains  separated  by  grain  boundaries  (a  microstructure)  does  not 
exist  in  glass.  Rather,  there  is  a  complete  lack  of  long-range  order,  but  short-  and 
intermediate-range  order  at  the  nano- structural  scale: 

•  Short-range  order:  Mostly  atom  to  atom  bond  lengths,  less  than  0.5  nm  and  bond  angles; 
characterized  by  radial  distribution  functions  (RDFs)  and  SAXS. 

•  Intermediate -range  order:  In  silica-based  glasses,  this  is  the  polymerization  of  the  silica  atomic 
tetrahedra  (one  Si  atom  surrounded  by  four  oxygen  atoms)  into  various  size  ring  structures  of 
joined  tetrahedra,  which  can  consist  of  4,  5,  6,  7,  8,  or  so  ring  groups  of  tetrahedra. 

Substitution  of  other  cations  (Na,  K,  Mg,  Ca,  etc.)  and  B  into  silica-based  glasses  can  have 
profound  effects  on  the  intermediate-range  order.  It  is  also  important  to  note  that  B  bonds  to 
three  oxygen  atoms  in  flat  triangles. 

•  Free-volume:  In  crystalline  materials,  using  the  theoretical  density,  it  is  straightforward  to 
calculate  the  atomically  unoccupied  free  space.  In  glasses,  this  unoccupied  atomic  space  is 
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referred  to  as  "free  volume,"  but  because  of  an  unknown  periodic  structure  it  is  extremely 
difficult  to  quantify  in  glasses.  The  free  volume  plays  a  critical  role  in  glass  densification  under 
stress/pressure. 

•  There  can  be  significant  complex  spatial  variations  of  defects,  free  volume,  atomic  structure, 
and  resulting  properties  at  the  nanoscale. 

•  It  is  the  current  wisdom  of  the  glass  community  that  the  short-  and  intermediate -range  order  at 
the  nanoscale  in  glasses  can  have  a  significant  influence  on  some  properties. 

4.1.3  Effect  of  Stress/Pressure 

Deformation  and  failure  in  ceramic  glasses  begins  with  reversible  to  irreversible  densification 

and/or  critical  cracks  nucleating  at  defects: 

•  High  pressure  can  have  significant  effect  on  coordination  of  Si  changing  from  typical  fourfold 
to  sixfold  coordination  of  oxygen  around  Si  -  (quartz-like  to  coesite-like  to  stishovite-like 
structures,  as  for  crystalline  quartz  shown  in  4). 

•  The  degree  of  polymerization  or  distribution  of  the  ring  structures  can  also  change  as  a  function 
of  pressure. 

•  The  common  consensus  is  that  the  ring  distribution  is  the  primary  control  of  some  properties; 
however,  at  very  high  pressures,  the  change  from  four  to  sixfold  Si  coordination  will 
significantly  influence  properties  as  well. 

•  Some  properties  of  fused  silica  (pure  Si02  glass)  are  anomalous,  e.g.,  a  negative  change  in 
shear  modulus  as  a  function  of  pressure.  Bando  et  al.  (3)  show  that  the  radius  of  curvature  of  a 
crack  in  glass  (figure  5)  can  be  about  1.5  nm,  suggesting  that  the  intermediate -range  order  can 
significantly  influence  crack  propagation. 
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Figure  5.  The  radius  of  curvature  of  a  crack  in  glass  after  Bando  (3). 
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4.1.4  Plasticity  in  Glass 


Ito  ( 4 )  has  presented  fairly  simple  methods  of  determining  the  brittleness  (figure  6)  or,  conversely, 
the  plasticity  of  glasses,  which  he  uses  to  suggest  that  the  brittleness  is  dependent  of  IRO  or  the 
distribution  of  ring  structures  seen  in  figure  7.  Table  2  lists  values  for  the  brittleness  parameters. 


Figure  6.  Brittleness  vs.  density  for  glasses  in  the  SiOo  and  /LLL-bascd 
glasses,  after  Ito  (4). 


(a)  (b) 

Figure  7.  (a)  Structure  of  soda  lime  glass  by  MD  simulation  where  the 

number  shown  is  the  ring  size  and  (b)  Number  of  network  rings 
in  soda  lime  (SL)  and  less  brittle  (LB)  glasses;  LB  glass  are 
more  polymerized  than  SL  glass  after  Ito  (4). 
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Table  2.  Brittleness  parameters. 


Glass 

Brittleness  Parameter  (/rm-1/2) 

Fused  Silica 

~10 

B2C>3  Based  glass 

~1 

Soda  Lime 

~5  -7 

The  brittleness,  therefore,  seems  to  be  dependent  on  the  deformation  and  fracture  behavior,  which 
depends  on  flow  and  densification  before  crack  initiation  and  on  the  bond  strength  of  the  network 
and  seems  to  decrease  with  a  decrease  in  density — a  free  volume  effect.  This  is  addressed  by 
Ito  (4)  in  the  same  paper. 

Conclusions  from  this  work  are  as  follows: 

•  Both  glasses  are  commercial  soda-lime  (SL)  based  glasses:  ( Na ,  K)20  —  (M g ,  Ca)0  -  Si02. 

•  The  less  brittle  (LB)  glass  appears  to  have  a  higher  polymerized  network  than  the  SL  glass,  i.e., 
a  significant  difference  in  the  ring  structure  distribution. 

•  IRO  at  the  nanoscale  seems  to  be  controlling  the  brittleness  of  these  glasses. 

4.2  Experimental  Results 

The  absorption/dissipation  of  deposited  energy  in  an  extreme  impact  event  depends  on  the  various 
deformation  and  failure/fracture  mechanisms  that  are  activated  during  the  event.  In  addition,  in  a 
multiscale  modeling  and  simulation  "Protection  Materials  by  Design"  approach  it  is  absolutely 
critical  to  experimentally  determine  the  key  properties  at  the  various  scales  to  validate  the 
theoretical  computational  results.  We  have  used  various  quasi-static  indentation,  edge-on-impact 
(EOI)  and  ballistic  impact  tests  for  this  purpose. 

4.2.1  Indentation 

Studies  on  the  deformation  and  fracture  of  glasses  and  AlON  using  a  spherical  500-/im-diameter 
diamond  indenter  was  recently  carried  out  by  Wilantewicz  (5).  Results  for  a  SL  (Starphire), 
boron  substituted  glass  (Borofloat),  fused  silica,  and  AlON  are  illustrated  in  figure  8  and  tables  3 
and  4,  there  are  significant  differences  in  the  deformation  and  fracture  behavior  of  these  materials. 
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Elastic  recovery  and  percent  elastic  recovery  as  a  function  of  the  maximum  indentation  load  for  AiON,  Starphire 
fused  silica  and  Borofioat glasses  using  a  500  mm  diameter  spherical  diamond  indenter,  ten  indentations  at 
each  load  were  made.  "Failure  Behavior  of  Glass  and  Aluminum  Oxynitride  {AION) 

Tiles  Under  Spherical  Ind enters”,  Trevor  E.  Wilantewicz.  ARL-TR45180  May  2010 


Figure  8.  Elastic  recovery  vs.  load  for  a  variety  of  glasses  after  Wilantewicz  (5). 


Table  3.  Elastic  recovery  in 

indentation  tests  at  200  N. 


Material 

%  Elastic  Recovery 

AION 

71 

Starphire 

73 

Borofioat 

79 

Fused  silica 

86 

Table  4.  Deformation  and  fracture  loads  after  Wilantewicz  (5). 


Onset 

All  Tests 

Onset  Ring 

All  Tests 

Onset  Radial 

All  Tests 

Dimpling 

Dimpled 

Cracking 

Ring  Cracked 

Cracking 

Radial  Cracked 

Material 

(N) 

(N) 

(N) 

(N) 

(N) 

(N) 

Staiphire  (tin) 

30 

30 

65 

75 

75 

100 

Starphire  (air) 

20 

30 

65 

100 

100 

125 

Borofioat  (tin) 

30 

35 

30 

45 

100 

200 

Silica  Glass 

75 

100 

20 

30 

65 

75 

AION 

35 

45 

45 

65 

40 

75 
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It  is  clear  that  these  materials  behave  in  significantly  different  ways.  The  onset  of  dimpling  is  the 
result  of  a  permanent  plastic  deformation.  The  normal  expectation  for  these  materials  is  that  as  a 
function  of  increasing  indentation  load  the  material  response  would  proceed  through  an  elastic 
regime,  then  through  a  plastic  regime  and  finally  into  a  cracking/fracturing  regime.  Table  5 
summarizes  the  observations.  Silica  glass  (fused  silica),  however,  behaves  in  a  dramatically 
different  way,  fracturing  prior  to  a  plastic  mechanism. 


Table  5.  Onset  of  elastic,  plastic,  and 
fracture  responses  after 
Wilantewicz  (5). 


Material 

Response 

Starphire: 

elastic  ->  plastic  — »  fracture 

Borofloat: 

elastic  — y  plastic  — >  fracture 

Silica  Glass: 

elastic  — y  fracture  — >  plastic 

AlON: 

elastic  — y  plastic  — >  fracture 

4.2.2  Edge-on-impact  Studies  of  Fused  Silica 

The  experimental  arrangement  is  illustrated  in  figure  9  and  a  series  of  time  resolved  photographs 
are  presented  in  figure  10. 


Plane 


Figure  9.  EOI  experimental  arrangement  after  Strassburger  et  al.  (6). 
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«  tl  #j  #1 


(a)  (b) 


Figure  10.  (a)  Illustration  of  a  series  of  EOI  tests  in  fused  silica  by  a  steel  solid 
cylinder  at  350  m/s  at  various  times;  first  and  third  rows  illustrate 
shaodwgraph  photos  in  plane  light  showing  damage;  second  and  fourth 
rows  are  photos  in  crossed  polarized  light,  which  show  the  propagation 
of  stress  via  a  photoelastic  effect,  and  (b)  illustrates  the  irregular  nature 
of  the  damage  front  due  to  the  presence  of  macro-defects  from  the 
same  test  after  Strassburger  et  al.  (6). 
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Table  6  lists  the  measured  velocities  of  the  longitudinal  waves,  transverse  waves  (shear  waves), 
and  crack  and  damage  front  velocities  in  fused  silica.  Note  that  the  longitudinal  wave  velocity 
for  fused  silica  is  5.93  km/s  and  the  shear  wave  velocity  is  3.77  km/s. 


Table  6.  Compilation  of  measured  wave,  crack,  and  damage  velocities  in 
fused  silica  after  Strassburger  et  al.  (6). 


Impact  velocity  (m/s) 

150 

260 
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Longitudinal  wave  speed  (m/s) 

shadowgraphs 

5975 

6076 

5823 

Longitudinal  wave  speed  (m/s) 

crossed  polarizers 

5814 

5796 

5491 

Transverse  wave  speed  (nr/s) 

shadowgraphs 

- 

3500 

3670 

Crack  velocity  (m/s) 

shadowgraphs 

2234 

2149 

2120 

Damage  velocity  (m/s) 

shadowgraphs 

5641 

5728 

5121 

4.2.3  Visualization  and  Analysis  of  Ballistic  Impact  Damage  and  Fragmentation  in 
Various  Glass  Plates 

In  this  series  of  experiments  Borofloat,  Starphire,  and  fused  silica  were  tested  in  a  standard 
ballistic  configuration.  The  plates  were  impacted  by  a  7.62-mm  AP  round,  and  a  solid  steel 
cylinder  inside  of  a  box  to  contain  all  of  the  resulting  fragments.  The  fragments  were  removed 
from  the  box  with  a  vacuum  and  then  sorted  by  sieves.  Figure  1 1  illustrates  the  experimental 
arrangement. 

Table  7  lists  the  details  of  the  various  tests  conducted  on  the  three  glass  types.  Note  that  the 
dimensions  of  the  Borofloat  glass  plate  used  in  test,  #  1 7742,  were  significantly  different  than  the 
others,  which  has  skewed  the  fragmentation  results  at  the  largest  sieve  size  2  mm.  Figure  12 
illustrates  a  series  of  very-high-speed  photographs  as  a  function  of  time  for  the  four  tests  listed  in 
table  7. 


Table  7.  Test  parameters  with  different  types  of  glass. 


EMI  Test  # 

Type 

Dimensions  (mm) 

Thickness  (mm) 

Projectile 

Impact  Velocity  (m/s) 

17742 

Borofloat 

149.4  x  149.7 

12.94 

cylinder 

1089 

17749 

Fused  silica 

101.65  x  101.67 

12.75 

cylinder 

1107 

17750 

Fused  silica 

101.62  x  101.62 

12.75 

7.62  AP 

824 

17751 

Starphire 

99.9  x  99.7 

10.06 

cylinder 

1115 

17 


(a)  (b) 


Figure  1 1 .  Schematic  of  (a)  ballistic  test  configuration,  and  (b)  target  after 
Strassburger  et  al.  (7). 


#  17751 :  soda  lime:  vp  = 
111 5 m/s:  Solid  cylinder 


#  17742:  Borofloat. 
solid  cylinder:  1089  m/s 


#  17749:fused  silica: 
solid  cylinder:  1100  m/s 


- 

* 


#  17750:fused  silica: 
7.62  AP:  800  m/s 


Figure  12.  Selection  of  high-speed  photographs  from  impact  on  various  glasses 
after  Strassburger  et  al.  (7). 
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The  observed  propagation  velocities  of  the  damage  zone  under  impact  of  the  steel  cylinder  are  all 
below  the  transverse  wave  velocities  of  the  materials,  as  seen  in  table  8,  which  compiles  the  wave 
and  fracture  velocity  data  determined  from  EOI  tests  along  with  data  from  the  actual  test  series. 
When  fused  silica  was  impacted  by  the  7.62-mm  AP  projectile,  the  formation  of  single  radial 
cracks  were  observed,  which  propagated  at  an  average  velocity  of  2394  m/s,  which  is  in  very 
good  agreement  with  the  crack  velocity  determined  for  fused  silica  from  EOI  tests. 


Table  8.  Compilation  of  wave  and  fracture  velocity  data  after  Strassburger  et  al.  (7). 


Glass  Type 

Longitudinal  Wave 
Velocity  (m/s) 

Transverse  Wave 

Velocity  (m/s) 

Terminal  Crack 

Velocity  (m/s) 

Damage 
Velocity  (m/s) 

Starphire 

5890 

3570 

1580 

3073® 

Borosilicate 

5543 

- 

2034 

2857® 

Fused  Silica 

6021 

3858 

2400 

3007®  2394t 

®  steel  cylinder,  v  =  1100  m/s 
1  AP  projectile,  v  =  824  m/s 


Figure  13  illustrates  the  fragmentation  data  of  the  various  glasses  listed  in  table  7.  The  data  in  the 
2-mm  range  should  be  ignored  as  the  plate  size  of  this  sample  was  much  different  than  the  others. 
The  fracture  and  damage  morphologies  using  the  solid  cylinders  compared  to  the  AP  round 
resulted  in  a  pronounced  wave-like  pattern  compared  to  Hertzian-like  radial  cracking  with 
differing  radial  patterns.  Analysis  of  the  fragmentation  characteristics  is  still  underway. 
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Figure  13.  Fragment  mass  distribution  from  sieve  analysis  for  tests  with  various 
glasses  after  Strassburger  et  al.  (7). 
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4.2.4  Nanoindentation  Studies  of  Fused  Silica 


Nanoindentation  experiments  were  conducted  at  ARL  on  fused  silica  specimens  using  an  MTS 
Nanoindenter  XP  operated  in  continuous  stiffness  measurement  (CSM)  mode.  A  spherical 
indenter  with  a  radius  of  3  /im  was  used  to  indent  to  depths  approaching  2  //in  at  a  constant  strain 
rate  of  0.05/s.  The  hardness  values,  H,  were  calculated  from  the  maximum  loads,  Pmax,  and  the 
contact  area,  A,  at  the  maximum  indentation  depth  where  H  =  The  elastic  modulus  values 

were  calculated  using  the  Oliver  and  Pharr  method  (16)  from  the  measured  unloading  stiffness,  S 
(which  is  equivalent  to  the  slope  of  the  initial  unloading  curve)  as  follows: 

S=^-EeffVA  (1) 

where  Eeff  is  a  function  of  the  Poisson’s  ratio  and  elastic  modulus  for  the  indenter  0/,,  Et)  and 
material  of  interest  (v,  E)  defined  as 


1  _  1  -  v2  1  -  v\ 

Eeff  E  +  Ei 


(2) 


Following  the  indentation  experiments,  the  residual  indents  were  examined  with  a  NanoSEM  600 
SEM  operated  in  low-vacuum  mode  (which  is  used  to  accommodate  non-conductive  specimens.) 

Figure  14  shows  the  typical  load-displacement  curves  for  maximum  displacements  ranging  from 
500  to  2000  nm.  Figure  15a  plots  the  measured  hardness  as  a  function  of  indentation  depth  and 
indicates  there  is  some  indentation  size  effect  over  the  range  considered.  The  standard  deviation 
is  largest  for  the  smallest  indentation  depths  (500  nm).  The  elastic  moduli  decrease  with 
increasing  maximum  indentation  depth  (figure  15b);  however,  the  error  in  measurement  does  not 
follow  a  trend  with  depth.  The  SEM  examination  gives  insight  into  the  indentation  size  effect. 
We  are  unable  to  resolve  a  residual  impression  for  the  specimens  indented  to  a  depth  limit  of  500 
nm,  which  indicates  the  response  is  mostly  elastic  at  small  depths.  Figures  16-18  show  SEM 
micrographs  for  specimens  indented  to  depths  of  1000,  1500,  and  2000  nm,  respectively.  There 
is  a  noticeable  residual  indent  in  figure  16;  however,  there  is  evidence  of  fracture.  At  greater 
indentation  depths,  radial  cracks  (figures  17a  and  18)  and  classic  cone  cracks  (figure  17b)  are 
visible.  Further  investigation  is  required,  but  the  lower  hardness  and  modulus  values  measured  at 
greater  indentation  depths  could  result  from  indentation  cracks.  In  some  brittle  material  systems, 
"pop-ins"  or  discrete  jumps  in  displacement  during  indentation  are  found  in  the  load-displacement 
curves.  However,  no  pop-ins  are  observed  in  this  series  of  nanoindentation  tests. 
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Figure  14.  Load  vs.  displacement  in  fused  silica. 
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Figure  15.  (a)  Hardness  and  (b)  modulus  variations  in  fused  silica. 
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Figure  16.  SEM  micrographs  of  fused  silica  indented  to  1000  nm. 


(a) 


(b) 


Figure  17.  SEM  micrographs  of  fused  silica  indented  to  1500  nm. 
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Figure  18.  SEM  micrographs  of  fused  silica  indented  to  2000  nm. 
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5.  Molecular  Dynamics  Modeling  of  Glass  Nanoindentation 


MD  methods  have  been  used  to  study  nanoindentation  for  a  number  of  material  systems,  i.e., 
metals  (17),  ceramics  (18),  glasses  (19),  and  energetic  materials  (20).  Length  scales  for  MD 
simulations  can  be  made  comparable  to  those  in  the  experiments  of  Nomura  et  al.  (19),  although 
the  time  scales  and  strain  rates  differ  significantly.  The  advantage  of  MD  simulations  over 
nanoindentation  experiments  is  the  capability  to  provide  atomistic  detail  of  numerous  properties, 
including  stress  distribution  and  structural  information. 

This  section  describes  preliminary  MD  simulations  of  the  nanoindentation  of  a  large-scale  fused 
silica  system.  All  the  simulations  were  performed  using  the  Large-scale  Atomic/Molecular 
Massively  Parallel  Simulator  (LAMMPS)  (21).  The  interatomic  potential  used  in  this  study  is  a 
pairwise  potential  recently  developed  using  force  matching  techniques  (11).  Force  matching  is  a 
fitting  approach  to  produce  classical  force  fields  from  trajectories  and  forces  obtained  from  ab 
initio  simulations  and  described  in  Izvekov  et  al.  (10).  The  development  of  such  a  potential  was 
necessary  due  to  inadequacies  of  existing  silica  potentials  in  describing  shock  properties,  as 
shown  in  figure  19  and  described  in  Gazonas  et  al.  (22),  which  compares  experimental  and 
predicted  shock  Hugoniot  information  using  existing  standard  classical  force  fields  and  that 
generated  using  force  matching.  Also,  unlike  the  standard  classical  force  fields,  the  force 
matched  potential  predicts  the  anomalous  densification  of  silica  at  high  pressures. 

The  amorphous  phase  of  Si02  (a -Si02)  was  generated  via  simulated  annealing.  In  general,  there 
are  two  means  of  annealing  a  system  from  a  melt  to  achieve  an  amorphous  phase:  by  decreasing 
the  temperature  in  discrete  intervals  or  by  continually  decreasing  the  temperature  during  the 
course  of  the  simulation.  In  both  cases,  the  system  is  simulated  in  the  constant  volume-constant 
temperature  (NVT)  ensemble,  with  periodic  boundaries  in  all  three  directions,  to  ensure  that  the 
final  density  matches  the  experimental  value.  For  this  study,  we  have  chosen  the  former  method, 
and  followed  the  annealing  schedule  described  in  Pedone  et  al.  (5).  The  temperature  of  the 
system  is  decreased  from  5000  K  in  500-K  increments,  resulting  in  a  cooling  rate  of  10  K/ps. 

The  system  in  this  study  consists  of  1,160,952  atoms  (386,984  a -Si02  molecules),  in  a  simulation 
cell  with  initial  dimensions  of  29.9  x  29.9  x  17.8  nm.  After  annealing  the  system,  indentation  is 
performed  by  a  spherical  indenter  with  a  radius  of  9  nm,  as  seen  in  figure  20.  The  indenter 
interacts  with  atoms  in  the  simulation  cell  via  a  force  of  magnitude 
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Figure  19.  Comparison  of  experimental  and  simulated  Flugoniot  curves  for  fused 
silica  (dL-Si02)'-  experiment  (•),  Morse  (8)  potential  (x).  BKS  (9) 
potential  (▲),  and  ARL  Force  Matching  method  (10,  11)  (M) . 
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Figure  20.  The  initial  configuration  of  the  fused  silica  system,  with  a 
representation  of  the  spherical  indenter. 
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(3) 


F  =  -K(r  -  R )2 

where  K  is  a  force  constant  (given  a  value  of  100.0  eV/A3  our  current  simulations),  R  is  the 
radius  of  the  indenter,  and  r  is  the  distance  from  the  atom  to  the  center  of  the  indenter.  The 
indenter  is  initially  centered  at  x  =  y  -  13.45  nm,  and  z  =  26.8  nm,  so  that  the  edge  of  the  indenter 
is  near  the  surface  of  the  fused  silica  system.  The  indenter  is  given  a  velocity  of  0.05  A/ps  (5 
m/s)  in  the  negative  2  direction.  The  system  can  be  unloaded  at  any  time  by  reversing  the 
velocity  of  the  indenter.  In  order  to  ensure  that  the  system  remains  stationary  during  indentation, 
approximately  33,000  atoms  (those  with  z  <  0.5  nm)  are  held  fixed.  Periodic  boundary 
conditions  are  implemented  in  the  x-  and  //-directions,  while  the  indented  surface  remains  free. 
The  simulation  is  performed  in  the  microcanonical  ensemble  (NVE),  with  a  timestep  of  2.0  fs. 

The  simulation  is  run  for  800,000  timesteps,  for  a  total  simulation  time  of  1.6  ns.  Figure  21 
shows  a  cross  section  of  the  atomic  configuration  at  four  different  times,  and  the  corresponding 
force  versus  distance  curve  is  shown  in  figure  22.  The  force  in  figure  22  is  the  ^-component  of 
the  force  on  the  indenter.  The  force  increases  with  distance  as  expected,  although  the  curve 
begins  to  plateau  near  the  end  of  the  simulation.  This  is  most  likely  due  to  the  proximity  of  the 
indenter  to  the  frozen  atoms  at  the  bottom  of  the  simulation  cell.  No  pileup  around  the  indenter  is 
observed,  and  no  cracking  occurs  during  loading.  The  system  is  unloaded  from  an  earlier  point 
in  the  simulation,  when  the  indenter  depth  is  approximately  4.0  nm,  to  avoid  artificial  effects  from 
the  frozen  atoms.  The  complete  loading  and  unloading  curve  is  shown  in  figure  23a.  The  cross 
sections  of  the  corresponding  atomic  configurations  at  the  maximum  load  and  after  complete 
unloading  are  shown  in  figures  23b  and  23c,  respectively.  From  these  images,  we  see  that  some 
elastic  recovery  has  taken  place.  The  hardness  is  estimated  as  the  maximum  load  divided  by  the 
projected  area  of  the  indent  and  is  determined  to  be  8.1  GPa.  This  is  slightly  less  than  the 
experimental  value  of  10  GPa  determined  by  Miyake  et  al.  (23);  however,  these  simulations  must 
be  repeated  on  larger  systems  to  ensure  that  the  size  of  the  sample  is  not  adversely  affecting  the 
hardness  results.  Additionally,  future  simulations  will  explore  the  effects  of  strain  rate  and 
indenter  size,  and  will  be  directly  compared  to  recent  experimental  results  described  in 
subsection  4.2.4  in  an  effort  to  determine  the  atomistic  mechanisms  driving  the  mechanical 
response  of  the  system  during  nanoindentation. 
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(a)  0.4  ns  (2  nm)  (b)  0.8  ns  (4  nm) 


Figure  21.  Cross  section  of  the  atomic  configuration  at  the  indicated  times.  The 
distances  in  parentheses  indicate  the  indentation  depth. 


Figure  22.  The  z-coniponcnt  of  the  force  on  the  indenter  plotted  against  the 
indentation  depth  corresponding  to  the  system  in  figure  21. 
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Figure  23.  (a)  Force  vs.  indenter  depth  for  a  complete  loading  and  unloading 

cycle.  Atomic  configuration  (b)  at  the  point  of  maximum  loading,  and 
(c)  after  complete  unloading. 


6.  Quantum  Mechanics  Modeling  of  Densification  and  Bulk  Modulus  of 
Silica  Versus  Pressure 


A  detailed  understanding  of  the  densification  process  and  structural  changes  in  amorphous  solids 
under  pressure  is  appealing  for  both  experimental  and  simulation  work.  One  of  the  interesting 
questions  is  about  the  nature  of  the  structural  transformation  between  low  and  high  density 
amorphous  phases.  To  model  the  structure  under  pressure  from  first  principles,  we  used  models 
with  different  density,  number  of  atoms,  and  different  ring  statistics.  We  used  two  different 
methods  to  generate  the  random  connected  networks  (1)  a  Monte-Carlo  bond  switching  model 
(72  and  1 14  atom  models)  and  (2)  MD  simulated  annealing  of  melted  silica.  The  ring  statistics 
could  be  described  by  plots  of  ring  size  distribution  seen  in  figure  24.  The  faster  the  quenching  of 
the  melt  is,  the  wider  is  the  distribution  of  the  ring  sizes.  Here  is  an  example  of  slow  quenching 
with  ring  distributions  from  4  to  8  member  rings,  which  corresponds  to  the  quartz-like  structure  of 
the  1 14  atom  model.  Two  other  models  with  72  and  192  atoms  have  wider  distribution  of  the  ring 
sizes  from  3  to  12  member  rings;  for  reference,  quartz  has  only  6  member  rings.  Three  to  four 
member  rings  have  a  low  concentration,  but  play  an  important  role  because  they  correspond  to  the 
most  reactive  sites.  The  angles  between  the  Si  —  O  —  Si  atom,  and  the  O  —  Si  —  O  atom 
distributions  also  convey  important  information  about  structural  changes  under  certain  pressure 
(figure  25). 
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Figure  24.  Ring  distribution  for  1 14  atom  model. 


Figure  25.  Two  types  of  angle  distribution  in  the  114  atom  model. 
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By  relaxing  the  internal  coordinates  under  compressive  or  tensile  pressures  at  normal  conditions, 
we  found  that  the  Si  —  O  —  Si  angle  ranges  between  130°  and  180°,  while  the  tetrahedral  units 
are  preserved.  This  result  is  in  good  agreement  with  previous  theoretical  and  experimental 
results  of  Mauri  et  al.  (24).  Clear  signs  of  structural  transformations  in  silica  under  the  pressure 
may  be  seen  from  the  calculated  mass  density  (figure  26).  We  did  the  calculations  of  density 
after  optimizing  the  shape  and  volume  of  the  unit  cell  without  projections  in  real  space.  The 
analysis  used  the  Projector  Augmented  Wave  (PAW)  method  implemented  in  the  Vienna  ab  initio 
simulation  package  (VASP)  code  described  in  Kresse  and  Hafner  (25). 
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Figure  26.  Density  of  fused  silica  as  a  function  of  pressure  for  72  atom  (Y),  1 14 
atom  (Y),  and  192  atom  (■jmodcls. 


In  all  three  models,  one  may  see  two  slopes  of  density,  which  might  be  related  with  different  types 
of  structural  changes.  One  occurs  up  to  20-30  GPa,  and  another  region  is  from  30-50  GPa. 

Some  signatures  of  the  two  phases  may  be  seen  from  the  x-ray  absorption  experiments  of  Sato 
and  Funamori  (72)  (figure  27). 
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Figure  27.  Experimental  density  of  fused  silica  after  Sato  and  Funamori  (72). 


To  understand  what  is  so  special  with  these  two  stages  of  structural  transformation  in  silica,  we 
did  angle  distribution  analysis  under  pressures  similar  to  that  depicted  in  figure  25.  It  turned  out 
that  in  the  first  region  there  is  not  much  change  in  O  —  Si  —  O  distribution,  but  discernible 
changes  in  Si  —  O  —  Si  distribution  indicating  that  up  to  20-30  GPa  there  is  a  squeezing  of  space 
between  tetrahedra  and  not  much  distortion  of  tetrahedrons.  Each  O  atom  has  two  nearest 
neighbors;  each  Si  atom  has  four  nearest  neighbors.  While  after  30-GPa  silica  becomes  very 
dense,  and  both  Si  —  O  —  Si  and  O  —  Si  —  O  angles  change,  revealing  transformation  in  both  the 
tetrahedra  and  the  space  between  them.  At  40-50  GPa  severe  distortion  of  the  tetrahedrons 
resulted  in  formation  of  sixfold-coordinated  Si  atoms  (figure  28).  The  DFT  observation  of  the 
sixfold-coordinated  atoms  confirms  the  assumption  suggested  in  Sato  and  Funamori  (72). 

A  manifestation  of  the  two  phases  of  silica  densification  might  exhibit  the  unusual  behavior  of  the 
elastic  constants,  particularly  the  bulk  modulus.  We  calculated  the  bulk  modulus  from  the 
stress-strain  relationship  for  1 14-atom  model.  The  bulk  modulus  of  fused  silica  generally 
increases  with  pressure,  but  unusual  behavior  of  the  pressure  dependence  up  to  20  GPa  may  be 
related  to  the  densification  of  the  space  between  the  tetrahedra  (figure  29).  A  second  region  after 
20  GPa  with  significant  increase  of  bulk  modulus  corresponds  to  the  densely  packed  tetrahedra. 
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Figure  28.  Fused  silica  structure  under  50-GPa  with  sixfold-coordinated  Si  atoms. 


Figure  29.  Bulk  modulus  of  fused  silica  as  a  function  of  pressure. 
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6.1  Force  Matching  Pair  Potentials  for  Borosilicate  Glasses 


To  generate  pair  potentials  for  borosilicate  glasses  we  used  the  force  matching  method  developed 
by  Izvekov  and  Rice  (11).  The  method  is  based  on  DFT  calculations  of  trajectories  in 
Bom-Oppenheimer  approximation  at  5000  K.  The  pair  potentials  for  S  —  B  and  B  —  O  generated 
using  the  method  are  shown  in  the  figure  30. 

Pair  Si  —  Si,  Si  —  O,  and  0  —  0  potentials  are  similar  to  those  used  for  pure  silica  MD 
calculations.  The  numerical  pair  potentials  include  to  a  certain  extent  many  body  interactions 
since  they  were  generated  based  on  force  matching  of  DFT  calculations. 


(a)  (b) 

Figure  30.  Pair  potentials  for  (a)  Si  —  B  and  (b)  O  —  B. 


7.  A  Perfectly  Matched  Layer  for  Peridynamics  in  Two  Dimensions 


In  this  section,  we  develop  a  peridynamic  method  for  modeling  nanoindentation  experiments  in 
fused  silica  and  other  chemically  substituted  glasses.  Originally  introduced  in  Silling  (26), 
peridynamics  is  a  non-local  formulation  of  elastodynamics,  which  can  more  easily  incorporate 
discontinuities  such  as  cracks  and  damage.  Derivatives  of  field  variables  in  the  classical 
continuum  model  are  replaced  by  integrals  over  a  small  neighborhood  of  microelastic  kernels, 
which  replaces  the  standard  constitutive  relation.  In  its  discretized  form,  an  elastic  solid  is 
treated  as  a  collection  of  particles  or  nodes,  each  connected  to  its  neighbors  by  breakable  bonds. 
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Bond  breakage  can  be  defined  to  occur  when  a  bond  is  stretched  past  some  predetermined  limit. 
The  end  result  is  a  method  capable  of  predicting  crack  growth  in  brittle  elastic  materials  (27-32). 
The  original  formulation  was  limited  to  materials  with  a  fixed  Poisson’s  ratio  of  0.25;  however, 
state -based  peridynamics  was  introduced  allowing  for  more  flexible  constitutive  relations  (33). 

As  will  become  clear  later,  state-based  peridynamics  allows  for  an  auxiliary  field  formulation, 
which  is  necessary  for  the  implementation  of  a  perfectly  matched  layer  (PML).  While  most 
peridynamics  work  has  focused  on  simulating  problems  with  free  or  fixed  boundary  conditions, 
there  are  applications  in  which  the  simulation  of  an  infinite  medium  may  be  useful,  such  as  crack 
propagation  in  a  halfspace  or  nanoindentation  problems.  Absorbing  boundary  conditions  are  one 
way  of  simulating  an  infinite  medium  as  any  impinging  waves  are  suppressed  so  they  do  not 
reflect  back  into  the  simulation.  A  PML  is  such  an  absorbing  boundary,  and  was  originally 
introduced  for  electromagnetic  simulations  (34,  35).  PMLs  differ  from  traditional  absorbing 
boundary  conditions  in  that  they  are  an  absorbing  layer  with  a  finite  width,  placed  between  the 
computation  region  of  interest  and  the  truncation  of  the  grid  or  mesh.  They  can  also  be  thought 
of  as  an  anisotropic  absorbing  material,  which  is  why  the  flexibility  of  a  state-based  peridynamics 
is  necessary.  A  PML  was  applied  to  one-dimensional  (1-D)  peridynamics  (36),  which  used  the 
results  of  Du  et  al.  (37)  to  formulate  an  auxiliary  field  equation.  This  approach  required  a  matrix 
representation  of  the  auxiliary  field,  which  may  be  memory  prohibitive  in  higher  dimensions. 

7.1  Two-dimensional  (2-D),  State-based  Peridynamics 

The  continuum  equation  of  motion  in  an  elastic  solid  can  be  stated  as 

Q 2  _ 

p— u  =  V  •  a  +  b  =  V  ■  (c  :  e)  +  b,  (4) 

where  p(x)  is  the  density,  u(x,  t)  is  the  displacement,  <x(x,  t)  is  the  stress  tensor,  e(x,  t)  is  the 
strain  tensor,  c  is  the  stiffness  matrix  for  plane  strain,  defined  by  Young’s  modulus  E,  Poisson’s 
ratio  v,  or  Lame  parameters  A  and  /i,  and  b(x)  is  a  body  force  (38).  (Throughout,  boldface  type 
denotes  a  vector  and  a  boldface  variable  with  an  overbar  denotes  a  tensor.)  Equation  4  is  a  local 
formulation  because  the  divergence  of  the  stress  (and  gradient  of  the  displacement  implied  in  its 
definition)  represents  a  local  operation  on  a  variable.  In  other  words,  the  action  of  V  ■  a  only 
depends  on  at  a  single  spatial  point.  In  problems  involving  discontinuities,  such  as  cracks,  the 
divergence  at  such  discontinuities  is  not  well  defined,  leading  to  numerical  implementation 
problems.  Peridynamics  proposes  replacing  V  •  a  with  a  nonlocal  operation  that  nonetheless 
also  represents  a  force.  Here,  we  use  the  state -based  peridynamics  (33),  rather  than  the  original 
bond-based  version  (26).  A  PML  application  requires  an  auxiliary  field  formulation,  as  it  is 
essentially  an  anisotropic  absorbing  material,  if  a  non-physical  one.  Consequently,  a  state-based 
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peridynamic  formulation  is  necessary  to  implement  the  required  constitutive  relations  in  the 
absorber.  The  state-based  peridynamics  uses  a  family  of  bonds  to  determine  a  given  force  rather 
than  a  single  bond  independently.  This  more  general  approach  allows  for  inelastic  behavior  and 
more  general  elastic  behavior,  and  is  governed  by 

f) 2  r  _  _ 

Pq^u=  (T[x,  t]  <x'  -  x)  -  T[x\  t]  (x  -  x'»  dV>  +  b)  (5) 

where  Hx  is  the  horizon  region,  defined  as  a  circle  centered  at  x  with  radius  5,  T[x,  t]  (x1  —  x)  is  a 
peridynamic  vector  state,  with  parameters  in  the  square  brackets  indicating  variables  that  act  as 
arguments  to  any  functions  that  define  the  vector  state  and  variables  in  the  angle  brackets  acting 
as  arguments  to  the  vector  state  itself.  In  the  state -based  formulation,  the  deformation  gradient, 
given  by 

F  =  I  +  uV,  (6) 

can  be  approximated  as  a  vector  state  as 


F[x,  t] 


cm(Y[x,t](t)®z)dvx, 


(7) 


where  C(|£|)  =  exp(—  |£|2  /S2)  is  a  shape  function,  taken  as  a  Gaussian  distribution  here  and 
with  the  horizon  T~LX  extended  so  that  the  shape  function  decays  to  an  arbitrary,  small  value,  taken 
here  as  10~6,  K  is  a  shape  tensor  given  by 


K[x,  t]  —  f  C(\£\)(Z®Z)dVx,,K 
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and  Y  is  a  deformation  vector  state  given  by 


(8) 


YM](0  =  v  + 


(9) 


with  ij  =  u[x',  t]  —  u[x,  t]  and  £  =  x'  —  x  (39). 

The  deformation  gradient  can  now  be  substituted  into  Hooke’s  law  and  strain-displacement 
relations,  giving  a  stress  term  cf  in  terms  of  u  in  plane  strain 

8 2 

Pfr2U  =  V  •  cr  =  V  •  (c  :  e) ,  (10) 

where 

e[x,  t\=l-  ( Vu  +  uV)  «  ^  (F[x,  t]  +  F[x,  t]T  -  21)  .  (1 1) 
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Ultimately,  the  peridynamic  vector  state  T  for  plane  strain  elasticity  is  given  by 


TM(£}  =  CflJDirMir'e 


(12) 


7.2  Auxiliary  Field  Formulation  and  PML  Application 

The  first  step  in  formulating  a  PML  is  to  construct  an  analytic  continuation  to  the  complex  plane, 
such  as  x  =  x  +  ig(x),  where  g(x)  is  a  given  function  describing  the  deformation  (40).  This 
mapping  has  the  effect  of  transforming  traveling  waves  of  the  form  elkx,  where  k  =  cu/c  is  the 
wave  number,  into  evanescent  waves  of  the  form  elkx e~k-kx\  thus  attenuating  such  waves  in  the 
PML  region.  Applying  a  PML  involves  substituting  for  spatial  derivatives  using 


d  Id 

dx  i  +  M  dx 

Ld 


(13) 


The  function  <fr(x)  defines  the  extent  of  the  PML  region,  i.e.,  when  o(x)  =  0  the  original  wave 
equation  is  obtained,  and  when  <p(x)  >  0,  traveling  waves  decay  exponentially.  Typically,  o{x) 
transitions  from  0  to  a  constant  value  using  a  smooth  function  to  prevent  numerical  reflections  at 
the  boundary  between  the  absorbing  and  computational  regions.  Before  applying  a  PML  directly 
to  the  peridynamic  equation,  equation  4  will  be  treated  so  that  the  PML  application  to 
peridynamics  will  be  clear.  It  is  convenient  to  convert  equation  4  to  the  Laplace  domain, 
assuming  e~st  time  dependence,  and  express  the  wave  equation  as  two  coupled  first  order  partial 
differential  equations,  the  first  in  u  and  the  second  in  s'ip  =  cf 


psu  =  V  ■  i/? 
sx/}  =  c  :  e, 


(14) 


where  the  Laplace  transform  of  a  variable  is  indicated  by  £{/}  =  /.  Expanding  equation  14  into 
components  gives  five  coupled  equations 
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(15) 


psux  = 

pSUy  = 

S^x  = 
S'lpy  = 
Slpr  = 


d  ~  d  ~ 

n  tyx  ~t~  rj 

ox  ay 
d  ~  d  ~ 
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d  d 

(A  +  2M)  -ux  +  \-u 

d  d 
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y 
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(  d  _  d  \ 
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We  can  expand  the  state -based  formulation  into  components  and  match  terms  to  equation  15. 
Following  this  approach  yields  a  viable  method  for  performing  PML  substitutions.  First,  the 
state -based  peridynamic  equations  listed  above  in  equations  5-12  can  be  written  explicitly  as 


psux[x,s]  =  C(|£|)  ('ij)x[x,s]k 


-i 


>  +  ( ^*[x',  s]kxx  +  W[x',  s]k  i 


+  /  cm  ('i/}x[x, s\kxy  +  ipT[x,  s]kyy  )^  +  (^X[x',s}kxy  +  ^T[x' ,  s]kyy  )  £y 
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psuy[x,  s]  =  C(|£|)  ( ^r[x,  s\kxx  +  ipy[x,  s}kyx  )  ix  +  ( Vv[x',  s\kxx  +  i[>y [x,  s\kyx  )  & 

Jhx 
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-l  i  J.  i  /"j.  rJ  7,-1  i  J.  7,-1 


s^x[x,s]  =  (A  +  2p) 


Yx[x,  s]ixkxx  +  Fx[x,  dVx,  -  1 


S'z/’y  [x, 
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„-i 


Fx[x,  +  yx[x,  s\£yk  d\ 4  -  1 
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+  A1  /  C-XICI)  ( y„[x,  +  yy[x,  s\£ykyx  )  dVx>, 
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(16) 


etc.  Though  no  derivatives  appear  in  equations  16,  the  correspondence  of  each  term  to  those  in 
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equation  15  is  apparent  and  the  PML  substitutions  can  be  made.  For  example,  the  first  equation 
in  16  can  be  rewritten  as 


(17) 


with  the  remaining  equations  following  similarly.  The  final  step  involves  simply  converting  back 
to  the  time  domain  and  implementing  a  forward  Euler  discretization  scheme  in  time,  and  the 
standard  one -point  integration  method  (32)  in  space. 

7.3  Results 

The  previous  method  was  numerically  implemented  using  a  forward  Euler  method  for  the 
temporal  discretization,  and  standard  one-point  integration  with  point-matching  in  space.  In  the 
following  example,  a  1  m-by-1  m  region  was  used  with  a  Gaussian  pulse  initial  condition  set  for 
the  x -directed  displacement.  For  the  PML,  a  Gaussian  ramp  was  used  as  a  smooth  transition  to  a 
constant  value,  set  at  50  s-1,  with  the  constant  region  having  a  width  0.1  m  and  the  ramp  a  width 
of  0.2  m.  Young’s  modulus  was  set  to  1  Pa,  Poisson’s  ratio  was  set  to  0.25,  and  the  density  was 
set  to  1  g/m2.  The  simulation  was  run  for  8  s.  Figure  31  shows  the  total  strain  energy,  summed 
over  the  entire  region,  as  a  function  of  time.  As  can  be  seen,  the  PML  absorbing  boundary  layer 
absorbs  impinging  waves  and  reduces  the  total  energy  by  five  orders  of  magnitude.  Figure  32 
shows  a  waterfall  plot  of  the  x-directed  displacement  along  the  y  =  0.5  m  line.  This  plot  also 
shows  the  decay  of  the  displacement  as  it  enters  the  PML  region. 


39 


Figure  3 1 .  Total  strain  energy  in  a  simulation  terminated  by  a  PML. 


Figure  32.  x-directed  displacement  at  y  =  0.5  m,  terminated  by  a  PML. 
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8.  Continuum  Equation-of-state  Model  Development 


Established  processes  to  develop  constitutive  models  for  use  in  large  scale  analysis  and  design 
simulations  are  heavily  influenced  by  the  vast  body  of  work  on  metals.  Typically,  a  model  is  split 
into  an  EOS,  relating  pressure,  density,  energy  and  temperature,  and  a  strength  model,  which 
relates  the  deviatoric  part  of  the  stress  tensor  to  the  time  history  of  the  strain  rate,  temperature, 
and  state  variables.  This  is  the  context  of  the  glass  EOS  in  Gazonas  et  al.  (22)  where  a 
polyamorphic  model  was  used,  in  conjunction  with  high  and  low  pressure  EOS  models,  to 
represent  the  permanent  volume  change  in  glass  at  high  pressures  (in  excess  of  10  GPa). 

The  model  assumptions  are  evaluated  by  comparison  of  simulation  results  with  velocity  profiles 
obtained  from  gas  gun  experiments  by  Alexander  et  al.  (13).  The  experimental  results  are  shown 
for  three  different  impact  velocities  in  figure  33a.  The  initial  velocity  ramp  from  0  to  0.2  km/s 
results  from  the  decreasing  modulus  with  pressure,  which  is  not  captured  by  the  current  model. 

At  velocities  of  approximately  0.6  km/s  there  is  a  kink  associated  with  the  polyamorphic 
transformation.  The  high  velocity  curve  (208  m/s  impact  velocity)  rises  abruptly  after  the 
transformation  point,  and  caps  at  a  velocity  consistent  with  the  initial  velocity  and  shock 
impedance  of  the  flyer.  The  transformation  is  nearly  over  driven  at  this  high  velocity,  in  that  there 
is  very  little  shift  between  the  two  vertical  sections  of  the  curve.  The  lower  curve  (137  m/s) 
shows  a  gradual  rise  associated  with  the  kinetics  of  the  transformation.  It  reaches  the  peak 
velocity  by  the  end  of  the  plot.  The  profile  from  the  intermediate  velocity  experiment  retains 
some  curvature  due  to  the  transformation  kinetics,  but  the  steady  velocity  is  attained  quickly. 

Results  from  using  the  initial  model,  with  separate  elastic -plastic  and  polyamorphic 
transformations,  in  gas  gun  simulations  yields  velocity  profiles  displaying  a  three-wave  structure 
(figure  33b  rather  than  the  two  waves  seen  in  the  experiments.  The  model  kinetics  were  assumed 
to  be  fast  in  these  simulations  to  accentuate  the  steps.  The  first  rise  in  the  profile  ends  with  the 
HEL  associated  with  the  elastic-plastic  transition.  The  second  rise  ends  with  the  volume  change 
from  the  polyamorphic  transformation,  and  the  third  rise  terminates  with  the  maximum  particle 
velocity.  From  this  comparison,  it  is  evident  that  the  elastic-plastic  and  polyamorphic 
transformations  must  occur  concurrently,  or  an  extra  step  will  appear  in  the  velocity  profile. 
Physically,  this  implies  that  the  atoms  move  to  accommodate  shear  strain  as  they  are  rearranging 
to  accommodate  the  volume  change. 
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Figure  33.  Velocity  profiles  from:  (a)  gas  gun  experiments  on  borosilicate  glass 
from  Alexander  et  al.  (13)  and  (b)  simulations  using  independent 
functions  for  the  elastic-plastic  transition  and  the  polyamorphic 
transformation. 


The  modeling  goal  is  to  couple  inelastic  deviatoric  deformation  and  volume  change.  There  is  no 
experimental  data  to  guide  the  functional  form  for  the  coupling,  so  a  simple  quadratic  flow 
potential  model  is  assumed: 


0=  0  =  \J  ao\  +p 2  —p(  A).  (18) 

Here,  ae  is  the  effective  stress,  p  is  the  pressure,  and  p  is  a  function  of  the  volume  fraction  of  the 
low  pressure  polyamorph,  A.  ct  is  a  parameter  regulating  the  relative  influence  of  the  deviatoric 
stress  compared  to  the  pressure. 

Following  concepts  used  in  metal  plasticity,  the  direction  of  inelastic  flow  is  assumed  to  be 
normal  to  the  flow  potential: 


d  inelas  A  ,  (19) 

OCT 

where  d ineias  is  the  inelastic  part  of  the  rate  of  deformation  tensor,  d,  which  is  decomposed, 
additively,  into  elastic  and  inelastic  parts  d  =  d,:/a.s  +  dineias.  The  elastic  part  is  also  a  function  of 
the  stress  tensor  through  the  elastic  moduli.  Together  these  provide  a  coupled  set  of  equations  for 
A  in  equation  19,  which  is  determined  by  an  iterative  procedure. 
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The  volume  fraction  of  the  low  density  polyamorph  is  assumed  to  evolve  through  a  kinetic 
relation 


A  =  /3(A  -  Aoc)3e(l/ 


(20) 


where  f3  and  v  are  parameters  and  Aoo  represents  the  equilibrium  fraction  of  the  low  density 
polyamorph.  The  equilibrium  phase  fraction 


(21) 


is  a  function  of  the  applied  loading,  p,  the  threshold  pressure  at  which  the  transformation  begins, 
pt ,  and  the  pressure  at  which  the  glass  is  fully  transformed  to  the  high  density  polyamorph,  pj . 

These  equations  were  implemented  into  a  finite  element  code  and  the  gas  gun  simulations  run. 
Figure  34a  shows  the  results  with  the  kinetic  parameter  set  high  (/ 3  =  100)  to  give  the  equilibrium 
response.  From  these  calculations  it  is  evident  that  a  two-wave  structure  is  produced,  indicating 
that  the  inelastic  deformation  and  polyamorph  transformation  are  occurring  simultaneously. 
Figure  34b  shows  the  results  with  a  kinetic  parameter  (/ 3  =  3)  set  to  reproduce  the  basic  features 
observed  in  the  experimental  data  plotted  in  figure  33a. 
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Figure  34.  Gas  gun  simulation  results  for  a  model  with  concurrent  elastic-plastic 


transition  and  polyamorphic  transformation  and  with  (a)  a  high  kinetic 
parameter  for  rapid  transformation  and  (b)  the  kinetic  parameter 
adjusted  to  resemble  experiments. 
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The  simulation  results  demonstrate  some  general  model  features  necessary  to  reproduce 
experimental  observations  from  high-pressure  gas  gun  experiments.  The  model  must  couple  the 
deviatoric  stress  and  the  pressure  for  the  polyamorphic  transformation,  and  a  kinetic  model  is 
necessary  to  introduce  time-dependent  rises  in  the  response.  Other  features  may  also  be 
necessary,  but  these  two  have  a  major  impact  and  provide  a  well  defined  starting  point  for  model 
improvement.  The  particular  functional  relations  assumed  for  this  model  are  conjecture  guided 
by  experience  in  modeling  metals.  Hence,  appropriate  data  are  needed  to  guide  model 
development.  While  it  will  be  possible  to  tune  a  kinetic  parameter  based  on  data  from  current 
experimental  methods,  existing  experimental  techniques  will  provide  little  information  on  the 
functional  form  for  the  kinetic  relation  or  the  flow  function  dependence  on  deviatoric  stress  and 
pressure.  The  alternative,  and  the  approach  pursued  in  the  program,  is  to  determine  the 
functional  relationships  and  parameters  through  multiscale  modeling  approaches. 


9.  A  Short-term  Conceptual  Project 


A  short-term  conceptual  project  to  determine  an  effective  experimental  and  theoretical  approach 
to  model  and  characterize  the  role  of  glassy  materials  in  resisting  ballistic  impact  was  conducted 
by  Richard  Lehman,  Professor  and  Chair  of  the  Department  of  Materials  Science  and 
Engineering,  Rutgers  University,  Piscataway,  NJ.  A  short  synopsis  of  the  report  was  delivered  to 
ARL  on  November  17,  2011  and  appears  below. 

A  number  of  features  of  the  glassy  state  are  thought  to  participate  in  the  ballistic  response  of  glass 
including: 

•  Structural  relaxation/viscoelasticity/strain  energy  equivalence  theory  (SEET) 

•  Short-  and  longer-range  atomic  structural  characteristics 

•  Free  volume  (compaction,  bulking) 

•  Cation  and  anion  coordination 

•  Bond  energies 

•  Other  nanoscale  order  characteristics  that  are  difficult  to  determine  experimentally  for  silicate 
glasses. 

A  number  of  expert  personnel  were  contacted  with  the  following  panel  discussion  outcomes: 

•  Structural  relaxation:  The  ability  of  glass  to  convert  to  a  liquid  with  little  structural  adjustment 
compared  to  crystalline  materials  may  enable  glass  to  respond  more  favorably  to  shaped  charge 
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assault  within  the  microsecond  time  periods  of  the  impact.  Viscoelastic  issues  will  be 
addressed. 

•  Modeling:  Structural  models  based  on  statistical  thermodynamics,  ring  structures,  molecular 
dynamics,  and,  anisotropic  finite  element  methods  (FEMs)  will  be  considered.  Experimental 
pair  distribution  functions  based  on  x-ray  data  will  be  discussed. 

•  Molecular  defects:  Glass  is  rich  in  molecular  defects  such  as  oxygen  hole  centers,  E,  bridging 
oxygen,  non-bridging  oxygen,  and  various  ring  defect  structure.  These  defects  may  play  an 
important  role  in  the  shaped  charge  behavior. 

•  Heterogeneous  structure:  Most  glasses  are  not  the  uniform  isotropic  material  at  the  mesoscale 
that  many  scientists  assume.  Phase  separation,  cation  clustering,  and  small-scale  density 
fluctuations  (as  evidenced  by  various  types  of  optical  scattering)  may  be  important  in  allowing 
glass  to  accommodate  high  strain  rates. 

•  Free- volume:  Glasses  contain  a  large  amount  of  free  volume  compared  to  their  crystalline 
counterparts. 

Summary  and  conclusions  of  the  study  include  the  following: 

•  MDs  can  be  expanded  into  the  mesoscale  region: 

-  Large  atom  arrays  (>109)  and  non-cubic  shapes,  e.g.,  high  aspect  ratio  cylinders. 

-  Use  surrogates  to  boost  scale. 

-  Time  scale  (microseconds)  cannot  be  collapsed  and  is  a  major  limit  on  modeling  time. 

•  FEMs  can  be  down-scaled  to  the  mesoscale  size  region: 

-  Requires  specialized  nonlinear  methods. 

-  Structural  elements  modeled  as  nonlinear  elements. 

-  Xi  Chen  at  Columbia  Univ.  has  relevant  experience. 

•  Medium-ranged  structural  data  obtained  with  scanning  tunneling  electron  microscopy  (STEM): 

-  Initial  material  characterization. 

-  Input  to  modeling  effort. 

-  Aid  in  interpretation  of  stress  wave  work. 

•  Stress  wave  characterization: 

-  Coherent  acoustic  phonon  spectroscopy  is  a  promising  dynamic  approach. 

-  The  time  scale  is  very  short. 

-  Experiments  can  be  adjusted  to  span  a  range  of  experimental  time  periods. 
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10.  Conclusions 


This  second-year  progress  report  on  multiscale  modeling  of  noncrystalline  ceramics  (glass)  has 
focused  on  establishing  the  framework  for  development  of  a  multiscale  computational 
methodology  for  optimizing  or  enhancing  the  performance  of  fused  silica  materials  not  yet 
synthesized.  A  more  immediate  research  objective  is  to  understand  why  certain  chemically 
substituted  fused  silica  materials  exhibit  enhanced  performance  in  the  defeat  of  SCJs  and  other 
ballistic  threats.  Conclusions  consistent  with  the  milestones  shown  in  the  five-year  roadmap 
shown  in  figure  3  are  as  follows: 

1.  Fused  silica  and  various  chemically  substituted  silica  specimens  have  been  delivered  to  ARL 
under  the  auspices  of  an  ongoing  cooperative  agreement;  these  specimens  have  been 
ballistically  tested,  and  will  serve  to  validate  future  computational  models  of  fused  silica. 

2.  The  structural  characteristics  of  silica  glasses  include  SRO  consisting  of  atom  to  atom  bond 
lengths  that  are  less  than  0.5  nm  and  bond  angles  characterized  by  RDFs  and  SAXS 
measurements.  IRO  in  silica-based  glasses  consists  of  the  polymerization  of  the  silica  atomic 
tetrahedra  (one  Si  atom  surrounded  by  four  0  atoms  into  ring  structures  of  joined  tetrahedra 
of  a  variety  of  sizes  ranging  from  4  to  8  groups  of  tetrahedra).  Chemical  substitution  of 

Na,  K,  Mg ,  Ca  and  other  atoms  can  have  a  profound  effect  on  IRO  and  ballistic  performance. 
Initial  free  volume  and  its  variation  with  pressure  control  densification  behavior,  but  this  is 
difficult  to  measure  or  quantify,  especially  at  elevated  pressure.  A  glass  brittleness  parameter 
discussed  by  Ito  ( 4 )  appears  dependent  on  the  deformation  and  fracture  behavior  which 
depends  on  flow  and  densification  before  crack  initiation  and  on  the  bond  strength  of  the 
network  and  seems  to  decrease  with  a  decrease  in  density. 

3.  A  series  of  experiments  on  Borofloat,  Starphire,  and  fused  silica  were  conducted  at 
Ernst-Mach  Institute,  which  showed  that  fracture  and  fragmentation  (figure  13)  of  these 
glasses  have  profoundly  different  macroscopic  response  (figure  12)  and  fracture  kinetics 
(table  8)  useful  for  computational  model  validation.  Solid  cylinder  impact  onto  fused  silica 
targets  results  in  a  greater  mass  of  finer  fragments  than  the  AP  round,  whereas  the  AP  round 
results  in  a  larger  mass  of  greater  than  2-mm-size  fragments  (figure  13). 

4.  A  series  of  nanoindentation  experiments  (indenter  radius  =  3  //m)  into  fused  silica  were 
conducted  at  ARL,  which  showed  that  both  hardness  and  modulus  decreased  with  increasing 
indentation  depth.  Also,  both  radial  and  cone  cracks  were  observed  in  indents,  which 
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exceeded  1500  nm  (figures  17  and  18).  The  nanoindentation  experiments  will  be  used  to 
validate  molecular  dynamics  and  peridynamics  models  of  fused  silica  and  other  chemically 
substituted  glasses. 

5.  Preliminary  MD  simulations  of  nanoindentation  (indenter  radius  =  9  nm)  into  fused  silica 
were  conducted  using  LAMMPS  (21)  using  a  pairwise  interatomic  potential  developed  using 
novel  force-matching  techniques  described  in  Izvekov  and  Rice  (11)  and  Izvekov  et  al.  (10). 
A  shock  Hugoniot  for  fused  silica  to  60  GPa  was  also  calculated  using  a  force-matched 
potential,  which  is  in  excellent  agreement  with  experimental  values  (figure  19). 

6.  First  principles  quantum  mechanical  methods  using  VASP  (25)  are  used  to  model 
densification  and  bulk  modulus  variations  with  pressure  in  fused  silica.  Ring  size 
distributions  (figure  24)  and  angle  distributions  for  Si  —  O  —  Si  and  O  —  Si  —  O  (figure  25) 
are  calculated,  and  future  work  will  consider  the  pressure  variation  of  these  distributions  and 
whether  they  can  be  validated  with  experimental  measurements  conducted  on  fused  silica 
under  extreme  pressure  using  diamond  anvil  cells.  Force-matched  potentials  were  also 
determined  for  a  borosilicate  glass  system  (figure  30). 

7.  A  2-D  peridynamic  model  was  developed  to  model  nanoindentation  and  dynamic  crack 
propagation  in  fused  silica  as  an  extension  of  the  recent  1-D  model  of  Wildman  and 
Gazonas  (36).  A  PML  was  added  to  permit  infinite  computational  domains  that  are 
encountered  in  solution  to  certain  boundary  value  problems. 

8.  Since  inelastic  deformation  and  a  polyamorphic  transformation  apparently  occurs 
simultaneously  in  some  glasses,  a  model  was  developed  to  account  for  this  coupled  behavior, 
which  captures  the  kinetics  of  the  polyamorphic  transformation  that  occur  in  plate  impact 
experiments  on  borosilicate  (compare  for  example  the  plate  impact  experiments  of  Alexander 
et  al.  (13)  in  figure  33a  with  simulation  results  in  figure  34b. 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


1- D 

2- D 

a-Si02 

Al 

AlON 

one-dimensional 

two-dimensional 

crystalline  quartz 

aluminum 

aluminum  oxynitride 

ARL 

a-Si02 

BKS 

B 

Ca 

CSM 

DAC 

DSI 

EOI 

EOS 

U.S.  Army  Research  Laboratory 
fused  silica  or  amorphous  quartz 

MD  potential  named  after  authors  (9) 

boron 

calcium 

continuous  stiffness  measurement 

diamond  anvil  cell 

Director’s  Strategic  Initiative 

edge-on-impact 

equation  of  state 

FEM 

finite  element  method 

FIB 

focused-ion-beam 

HEL 

IRO 

LAMMPS 

Hugoniot  elastic  limit 
intermediate-range  order 

large-scale  atomic/molecular  massively  parallel  simulator 

LB 

less  brittle 

MD 

molecular  dynamics 

Na 

sodium 

MEDE 

materials  in  extreme  dynamic  environments 

NVT 

0 

PAW 

constant  volume-constant  temperature 

oxygen 

projector  augmented  wave 

PML 

perfectly  matched  layer 

RDF 

radial  distribution  function 

SAXS 

SCJ 

SEET 

SEM 

small-angle  x-ray  scattering 
shaped-charge  jet 
strain  energy  equivalence  theory 
scanning  electron  microscopy 
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List  of  Symbols,  Abbreviations,  and  Acronyms  (Continued) 


Si 

silicon 

SiC 

silicon  carbide 

SL 

soda  lime 

SRO 

short-range  order 

STEM 

scanning  tunneling  electron  microscopy 

VASP 

Vienna  ab  initio  simulation  package 

WMRD 

Weapons  and  Materials  Research  Directorate 
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1  JOHNS  HOPKINS  UNIV 
DEPT  OF  MECH  ENGRG 
K  T  RAMESH 
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4  RDRLD 

C  CHABALOWSKI 
J  CHANG 
O  OCHOA 
R  SKAGGS 
BLDG  205 

2800  POWDER  MILL  RD 
ADELPHI  MD  20783-1197 


ABERDEEN  PROVING  GROUND 

86  DIR  USARL 
RDRL  WM 
B  FORCH 
S  KARNA 
J  MCCAULEY 
P  PLOSTINS 
P BAKER 
RDRL  WML 
J  NEWILL 
M  ZOLTOSKI 
RDRL  WML  B 
I BATYREV 
S  IZVYEKOV 
BRICE 

R  PESCE  RODRIGUEZ 
D  TAYLOR 
N  TRIVEDI 
N  WEINGARTEN 
RDRL  WML  D 
P  CONROY 
M  NUSCA 
RDRL  WML  E 
P  WEINACHT 
RDRL  WML  F 
D  LYON 
RDRL  WML  G 
M  BERMAN 
W  DRYSDALE 
RDRL  WML  H 
D  SCHEFFLER 
S  SCHRAML 
B  SCHUSTER 
RDRL  WMM 
J  BEATTY 
R  DOWDING 
J  ZAB INSKI 
RDRL  WMM  A 
J  TZENG 
E  WETZEL 
RDRL  WMM  B 


T  BOGETTI 
B  CHEESEMAN 
C  FOUNTZOULAS 
G  GAZONAS 
D  HOPKINS 
R  KARKKAINEN 
BLOVE 
PMOY 
B  POWERS 
C  RANDOW 
TSANO 
F  TAVAZZA 
M  VANLANDINGHAM 
R  WILDMAN 
C  YEN 

RDRL  WMM  C 
J  LA  SCALA 
RDRL  WMM  D 
ECHIN 
KCHO 

RDRL  WMM  E 
J  ADAMS 
M  COLE 
T  JESSEN 
J  LASALVIA 
P  PATEL 
J  SANDS 
J  SINGH 
RDRL  WMM  F 
L  KECSKES 
H  MAUPIN 
RDRL  WML  G 
J  ANDZELM 
A  RAWLETT 
RDRL  WMP 
B  RINGERS 
S  SCHOENFELD 
RDRL  WMP  B 
M  GREENFIELD 
C  HOPPEL 
R KRAFT 
M  SCHEIDLER 
T  WEERASOORIYA 
RDRL  WMP  C 
R  BECKER 
S  BILYK 
T  BJERKE 
D  CASEM 
I  CLAYTON 
B  LEAVY 
M  RAFTENBERG 
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S  SATAPATHY 
S  SEGLETES 
RDRL  WMP  D 
R  DONEY 
D  KLEPONIS 
I RUNYEON 
B  SCOTT 
H  MEYER 
RDRL  WMP  E 
M  BURKINS 
RDRL  WMP  F 
A  FRYDMAN 
N  GNIAZDOWSKI 
R  GUPTA 
RDRL  WMP  G 
N  ELDREDGE 
D  KOOKER 
S  KUKUCK 
G  R  PEHRSON 
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